1. Field of the Invention
Equation Chapter 1 Section 1
This invention relates generally to the field of geophysical prospecting. More particularly, the invention relates to the field of seismic data processing. Specifically, the invention is a method for seismic wavefield extrapolation using variable extrapolation step size and phase-shifted linear interpolation.
2. Description of the Related Art
The use of three-dimensional (3D) seismic methods has resulted in increased drilling success in the oil and gas industry. However, 3D seismic methods are computationally expensive. A crucial point in processing 3D seismic data is the migration step, both because of its 3D nature and the computational cost involved. Seismic migration is the process of constructing the reflector surfaces defining the subterranean earth layers of interest from the recorded seismic data. Thus, the process of migration provides an image of the earth in depth or time. It is intended to account for both positioning (kinematic) and amplitude (dynamic) effects associated with the transmission and reflection of seismic energy from seismic sources to seismic receivers. Although vertical variations in velocity are the most common, lateral variations are also encountered. These velocity variations arise from several causes, including differential compaction of the earth layers, uplift on the flanks of salt diapers, and variation in depositional dynamics between proximal (shaly) and distal (sandy to carbonaceous) shelf locations.
Seismic depth migration has traditionally been performed through the application of Kirchhoff methods. However, because of recent advances in computing hardware and improvements in the design of efficient extrapolators, methods that are based on solutions of the one-way wave equation are becoming increasingly popular. Downward wave extrapolation results in a wavefield that is an approximation to the one that would have been recorded if both sources and receivers had been located at a depth z. Downward extrapolation of the seismic wavefield is a key element of wave equation based migrations because it determines the image quality and computational cost. Extrapolation by finite difference methods is utilized for its ability to handle lateral velocity variations. It can be classified into explicit and implicit methods.
A number of the algorithms used in two-dimensional (2D) depth migration have not been successful in the 3D case. For example, implicit finite-difference extrapolation methods have the advantage of being unconditionally stable, but the disadvantage of being difficult to extend to the 3D case. The most common 3D implicit method is based on the operator splitting into alternating direction components. The splitting errors that these methods exhibit in the 3D case are translated into non-circularly symmetric impulse responses, which become unacceptable for dips higher than 45°. The splitting methods have errors that depend significantly on reflector dip and azimuth and thus have problems with reflector positioning errors. Similarly, two-pass methods have problems handling lateral velocity variations. In contrast, explicit extrapolation methods that approximate the extrapolation operator as a finite-length spatial filter are easily extended to the 3D case. The difficulty with explicit extrapolation is that the stability condition is not automatically fulfilled. The stability condition is that no amplitude at any frequency will grow exponentially with depth. Stability must be guaranteed by careful design of the extrapolation operators.
Three-dimensional seismic wavefields may be extrapolated in depth, one frequency at a time, by two-dimensional convolution with a circularly symmetric, frequency- and velocity-dependent operator. This depth extrapolation, performed for each frequency independently, lies at the heart of 3D finite difference depth migration. The computational efficiency of 3D depth migration depends directly on the efficiency of this depth extrapolation. For these techniques to yield reliable and interpretable results, the underlying wavefield extrapolation must propagate the waves through inhomogeneous media with a minimum of numerically induced distortion over a range of frequencies and angles of propagation.
Holberg, 0., “Towards optimum one-way wave propagation”, Geophysical Prospecting, Vol. 36, 1988, pp. 99–114, discloses the first explicit depth extrapolation with optimized operators. Holberg (1988) proposes a technique, solely for 2D depth migration, by generalizing conventional finite-difference expressions in the frequency-space domain. This technique yields optimized spatial symmetric convolutional operators, whose coefficients can be pre-computed before migration and made accessible in tables. The ratio between the temporal frequency and the local velocity is used to determine the correct operator at each grid point during the downward continuation, by fitting their spatial frequency response to the desired phase shift response over a range of frequencies and angles of propagation to control numerical distortion. The Holberg (1988) technique can be made to handle lateral velocity variations, but only applies to 2D migration.
Blacquiere, G., Debeye, H. W. J, Wapenaar, C. P. A., and Berkhout, A. J, “3-D table driven migration”, Geophysical Prospecting, Vol. 37, 1989, pp. 925–958, extends the method disclosed in Holberg (1988) to 3D migration. The wavefield extrapolation is performed in the space-frequency domain as a space-dependent spatial convolution with recursive Kirchhoff extrapolation operators based on phase shift operators. The optimized operators are pre-calculated and stored in a table for a range of wavenumbers. The extrapolation is performed recursively in the space domain, so both vertical and lateral velocity variations can be handled. The Blacquiere et al. (1989) method is accurate, but is a full 3D method and hence computationally expensive.
Hale, D., “3-D depth migration via McClellan transformations”, Geophysics, Vol. 56, No. 11 (November, 1991), pp. 1778–1785, introduces a more efficient 3D scheme based on the McClellan transform, which gives numerically isotropic extrapolation operators. This is commonly referred to as the Hale-McClellan scheme. Given the coefficients of one-dimensional frequency- and velocity-dependent filters similar to those used to accomplish 2D depth migration, McClellan transformations lead to an algorithm for 3D depth migration. Because the coefficients of two-dimensional depth extrapolation filters are never explicitly computed or stored, only the coefficients of the corresponding one-dimensional filters are required, leading to a computationally more efficient scheme. However, the Hale (1991) method has numerical anisotropy.
Soubaras, R, “Explicit 3-D migration using equiripple polynomial expansion and Laplacian synthesis”, 62nd Ann., Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1992, pp. 905–908, improves the Hale-McClellan scheme. Soubaras (1992) uses an expansion in second-order differential operators instead of the McClellan transform. The Soubaras method also uses the Remez algorithm to design the coefficients for both the extrapolation operator and for the differential operators. This method avoids the numerical anisotropy to a large degree and is comparable in computational cost to the Hale-McClellan scheme. The Soubaras approach takes advantage of operator circular symmetry, as do the McClellan transformations, but avoids the computation of a 2D filter approximating the cosine of the wavenumbers. The approach defines a Laplacian operator, approximates the Laplacian by the sum of two 1D filters approximating the second derivatives, and approximates the exact extrapolation operator by a polynomial. The second derivative operators and the polynomial expansion are both calculated by the Remez exchange algorithm.
Sollid, A., and Arntsen, B., “Cost effective 3D one-pass depth migration”, Geophysical Prospecting, Vol. 42, 1994, pp. 755–716, makes the Soubaras (1992) scheme more cost effective. Sollid and Arntsen (1994) used frequency adaptive optimized derivative operators. The expansion in second order differential operators is used, but the expansion coefficients and differential operators are designed using a least squares approach rather than the Remez algorithm. A set of variable length second order differential operators for each wavenumber has different spectra and lengths to ensure that the resulting wave extrapolation is both accurate and efficient.
Biondi, B., and Palacharla, G., 1995, 3D depth migration by rotated McClellan filter: Geophysical Prospecting, 43, 1005–1020, improves the accuracy of the Hale-McClellan scheme. Biondi and Paracharla (1995) rotate the 9-point McClellan filter by 45° and then average the rotated filter with either of the original 9-point or 17-point McClellan filters. The averaged filter images steep dips without dispersion of the wavefield and without much additional computational cost.
Mittet, R., 2002, “Explicit 3D depth migration with a constrained operator”, 72nd Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, pp. 1148–1151, discloses a constrained explicit operator method which is a modification of Blacquiere et al. (1989) to make the scheme more computationally efficient. The number of independent operator coefficients is constrained to reduce the number of computer floating point operations required, thus increasing computer efficiency. The innermost coefficients in the core area of the extrapolation operator are computed in a standard fashion. The remaining outermost coefficients in the operator, related to very high dip and evanescent wave propagation, change only as a function of radius and are constant within radial intervals.
All of the methods discussed above are, however, still computationally expensive. Thus, a need exists for an explicit depth extrapolation method for 3D seismic migration that is more computationally efficient.